Libraries:
# libraries -------------------
library(ggplot2)
library(vegan)
library(ggvegan)
library(tidyverse)
Set working directory
knitr::opts_knit$set(root.dir = '/Users/user/Desktop/Data/Regen_RProj/')
Functions:
# functions -------------------
source("Functions/veganCovEllipse.R")
mali_wider_log <- as.data.frame(lapply(mali_wider, as.logical))
ma.li <- as.data.frame(colSums(mali_wider_log)) %>%
rownames_to_column(var = "Species")
# Species from the other 5% calcs
# In NH/ME - no LYQU, CEOR, PAQU, CEAM, ARME, LYBO, RHGL
# For this data set , 5% is ~11 plots:
# Would add: SMGL
# If I stuck to at least 11, would lose: ACRU, COPE, PRSE, PIST, POUN11
rm(mali_wider_log, ma.li)
# 3 dimensions
mali_3d <- metaMDS(rel_mali5,
distance = 'bray',
k = 3,
trymax = 20,
autotransform = FALSE)
## Run 0 stress 0.1395122
## Run 1 stress 0.1465263
## Run 2 stress 0.1369332
## ... New best solution
## ... Procrustes: rmse 0.1183403 max resid 0.3726543
## Run 3 stress 0.1391444
## Run 4 stress 0.1458937
## Run 5 stress 0.1530509
## Run 6 stress 0.1358958
## ... New best solution
## ... Procrustes: rmse 0.05522124 max resid 0.1657867
## Run 7 stress 0.1369326
## Run 8 stress 0.1358959
## ... Procrustes: rmse 0.0005646437 max resid 0.001453757
## ... Similar to previous best
## Run 9 stress 0.1465264
## Run 10 stress 0.136823
## Run 11 stress 0.1358959
## ... Procrustes: rmse 0.0002755069 max resid 0.0008064717
## ... Similar to previous best
## Run 12 stress 0.1341827
## ... New best solution
## ... Procrustes: rmse 0.0594534 max resid 0.2139219
## Run 13 stress 0.1376197
## Run 14 stress 0.1376205
## Run 15 stress 0.1391435
## Run 16 stress 0.1395122
## Run 17 stress 0.135896
## Run 18 stress 0.1458936
## Run 19 stress 0.1369321
## Run 20 stress 0.1405414
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 20: stress ratio > sratmax
# converges, stress of 0.13 - 0.15
# 2 dimensions
mali_2d <- metaMDS(rel_mali5,
distance = 'bray',
k = 2,
try = 40,
trymax = 40,
autotransform = FALSE)
## Run 0 stress 0.2239971
## Run 1 stress 0.2187381
## ... New best solution
## ... Procrustes: rmse 0.1215493 max resid 0.3225023
## Run 2 stress 0.224739
## Run 3 stress 0.2047202
## ... New best solution
## ... Procrustes: rmse 0.1516674 max resid 0.3370896
## Run 4 stress 0.2066462
## Run 5 stress 0.2360426
## Run 6 stress 0.2064071
## Run 7 stress 0.2216095
## Run 8 stress 0.2097745
## Run 9 stress 0.2274517
## Run 10 stress 0.2191102
## Run 11 stress 0.2047202
## ... New best solution
## ... Procrustes: rmse 0.00158856 max resid 0.005415167
## ... Similar to previous best
## Run 12 stress 0.2150753
## Run 13 stress 0.2080347
## Run 14 stress 0.2211264
## Run 15 stress 0.3659444
## Run 16 stress 0.2150752
## Run 17 stress 0.2289216
## Run 18 stress 0.2272952
## Run 19 stress 0.2245979
## Run 20 stress 0.2234159
## Run 21 stress 0.2132488
## Run 22 stress 0.224739
## Run 23 stress 0.2150117
## Run 24 stress 0.2247403
## Run 25 stress 0.2232724
## Run 26 stress 0.2257935
## Run 27 stress 0.2203333
## Run 28 stress 0.2097744
## Run 29 stress 0.2132487
## Run 30 stress 0.2297719
## Run 31 stress 0.2185977
## Run 32 stress 0.2237861
## Run 33 stress 0.2301115
## Run 34 stress 0.2261155
## Run 35 stress 0.2308032
## Run 36 stress 0.2185277
## Run 37 stress 0.2150117
## Run 38 stress 0.2180135
## Run 39 stress 0.2132487
## Run 40 stress 0.2146601
## *** Best solution repeated 1 times
# stress is too high to represent in 2 dimensions
ma.veg_enr <- envfit(mali_3d, ma.meta.data_veg)
ma.veg_spc.envr <- envfit(mali_3d, rel_mali5)
# site data -------------------
ma.site.scrs <- as.data.frame(scores(mali_3d, display = "sites"))
ma.site.scrs <- cbind(ma.site.scrs, Treat_Type = ma.meta.data_veg$Treat_Type)
ma.site.scrs <- cbind(ma.site.scrs, Region = ma.meta.data_veg$Region)
# species data -------------------
ma.spp.scrs <- as.data.frame(scores(ma.veg_spc.envr, display = "vectors"))
ma.spp.scrs <- cbind(ma.spp.scrs, Species = rownames(ma.spp.scrs))
ma.spp.scrs <- cbind(ma.spp.scrs, pval = ma.veg_spc.envr$vectors$pvals)
sig.ma_spp.scrs <- subset(ma.spp.scrs, pval<=0.05)
# envr data -------------------
ma.envr.scrs <- as.data.frame(scores(ma.veg_enr, display = "vectors"))
ma.envr.scrs <- cbind(ma.envr.scrs, env.variables = rownames(ma.envr.scrs))
ma.envr.scrs <- cbind(ma.envr.scrs, pval = ma.veg_enr$vectors$pvals)
ma.sig_envr.scrs <- subset(ma.envr.scrs, pval<=0.05)
library(plotly)
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
fort.mali3 <- fortify(mali_3d) %>%
filter(score == "sites")
plot_ly(x = fort.mali3$NMDS1, y = fort.mali3$NMDS2, z= fort.mali3$NMDS3, type = "scatter3d", mode = "markers", color = ma.meta.data_veg$Treat_Type, symbols = ma.meta.data_veg$Region)
summary(as.factor(ma.meta.data_veg$Treat_Type))
## Control FallRx Harvest MowRx SpringRx
## 5 3 6 3 6
ma.tt <- adonis2(rel_mali5 ~ Treat_Type*Region, data = ma.meta.data_veg)
print(ma.tt)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = rel_mali5 ~ Treat_Type * Region, data = ma.meta.data_veg)
## Df SumOfSqs R2 F Pr(>F)
## Treat_Type 4 1.7888 0.26616 1.8764 0.002 **
## Region 1 0.7921 0.11786 3.3236 0.001 ***
## Treat_Type:Region 2 0.5649 0.08406 1.1852 0.261
## Residual 15 3.5749 0.53192
## Total 22 6.7207 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#each are significant, but the interaction is not